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We present large-scale Molecular Dynamics simulation results for granular gases of particles with 
the velocity-dependent restitution coefficient, as it follows from the simplest first-principle model of 
viscoelastic spheres. We analyze both cases of force-free gases and homogeneously heated gases. We 
develop a theory for the velocity distribution for heated gases. We explored evolution of temperature 
and the Sonine coefficient, which characterize the form of the velocity distribution function. We 
observed that for not large dissipation the simulation results are in excellent agreement with both, the 
existing theory for gases in a homogeneous cooling state and with the novel theory for heated gases. 
At the same time a noticeable discrepancy between theory and MD results have been detected for 
large dissipations. We analyze the adequacy of the simplified model of a constant effective restitution 
coefficient for uniformly heated gases of viscoelastic particles. We conclude that while this model 
may be used to describe temperature of a gas, it fails to describe quantitatively both temperature 
and the velocity distribution function. 

PACS numbers: 81.05.Rm, 05.20.Dd, 05.40.2a 



I. INTRODUCTION 

Granular materials are abundant in nature and exam- 
ples range from sands and powders on Earth EtQ to 
more dilute systems, termed as granular gases [aTMl '^^ 
space; astrophysical objects, like planetary rings, may 
be mentioned in this respect [l^l • These systems are 
of a great scientific and industrial importance. Depend- 
ing on the average kinetic energy of individual grains, 
granular materials can be in a solid, liquid or gaseous 
state. In this paper, we address a low-density granu- 
lar gas - a system comprised of macroscopically large 
number of macroscopic grains that move ballistically be- 
tween inelastic pair-wise collisions. The inelastic nature 
of inter-particle collisions is the main feature which dis- 
tinguishes a granular gas from a molecular gas and gives 
rise to many astonishing properties of dissipative gases. 
If no external forces act on the system it evolves freely 
and permanently cools down. During the first stage of 
a granular gas evolution its density remains uniform and 
no macroscopic fluxes present. This state, which is called 
a homogeneous cooling state (HCS) |Tl|, is however, 
unstable with respect to fluctuations and in later times 
clusters and vertexes develop. The system then continue 
to loose its energy in the inhomogeneous cooling state 
(ICS) [T3 - [l9j . If energy is injected into a granular gas to 
compensate the loss of kinetic ener gy, the system settles 
into a nonequilibrium steady-state (20l . [2]| . 

The dissipation of energy in pair-wise collisions is 
quantified by the restitution coefficient e. 
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where v^'j = V2 — v{ and V12 = V2 — vi are the relative 
velocities of two particles after and before a collision, and 



e is a unit vector connecting their centers at the collision 
instant. The post-collision velocities are then related to 
the pre-collision velocities Vi and V2 as follows, e.g. Q: 
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To date, most studies of the HCS and ICS have fo- 
cused on the case of a constant restitution coefficient [iTI - 
fioj . This ass ump tion, however, contradicts experimental 
observations j22 - l24l |. along with basic mechanical laws 
[25I [26j , which show that e does depend on the impact 
velocity [23, [26l - [29| . This dependence may be obtained 
solving equation of motion for colliding particles with the 
explicit account for the dissipative forces acting between 
the grains. The simplest first-principle model of inelas- 
tic collisions assumes viscoelastic properties of particles 
material, which result in viscoelastic inter-particle force 
[2?! and finally in the restitution coefficient [1^, [2^ : 
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the numerical coefficients are Ci 
The elastic constant 
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is a function of the Young's modulus Y and Poisson ra- 
tio 1^ and mass m; the constant A quantifies the viscous 
properties of the particles' material. While the last equa- 
tion is valid for not very large inelasticities, in the recent 
work of Schwager and Poschel [s^l an expression for e, 
which is valid for highly inelasticity and accounts for the 
delayed recovery of the shape of colliding particles has 
been derived: 
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Here, hi = 0, h2 = — Ci, /13 = 0, = C2, and the other 
numerical coefficients up to iV/j = 20 are given in Ref. 

In the latter equation u{t) = T{t)/To {To = T(0)) 
is the reduced temperature, with the usual definition of 
temperature 
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Here n is the number density of the gas, v is a particle ve- 
locity, /(v,t) - the time-dependent velocity distribution 
function and C12 is the dimensionless relative velocity, 
defined as C12 = V12/VT, where vrit) = y^2T{t)/m is 
the thermal velocity. The dissipative constant 5 reads. 
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so that the first few terms in Eq. ([5]) are identical to 
those in Eq. ([3]). Note that i5 = corresponds to an 
elastic system, e = 1 and the dissipation in the system 
increases with increasing S. Furthermore, e — > 1 when the 
reduced temperature u{t) decreases, that is, the effective 
restitution coefficient increases as a system cools down. 
The dependence of e on the dissipation parameter S is 
shown in Fig. [T] for u{t) = 1 and the average value of 
the reduced normal velocity (|ci2 • e|) = 4/\/2tt. As it 
may be seen from the figure, the restitution coefficient, 
given by Eq. ([5]) demonstrates a good convergence for 
d < 0.5. In what follows we use the number of terms in 
the expansion ([5]) equal to Nk = 20, which guarantees 
an accurate description of e up for S < 0.5, that is, for 
rather large dissipation. 

The velocity distribution function in granular gases 
usually deviates from the Maxwellian distribution [l^, 
figl . IsimSj and may be represented in terms of a Sonine- 
polynomial expansion Q : 
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where the Sonine polynomials read for d = 3: 

(-l)'=(l/2-t-p)! 
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Hence evolution of ./(v,t) is completely determined by 
the time dependence of the Sonine coefficients ap(t). The 
first Sonine coefficient is trivial, oi = 0, as it follows from 
the definition of temperature, e.g. [1| while the elastic 
Maxwellian limit corresponds to Cp = for all p > 1. 

It has been shown recently that for granular gases with 
a constant e the expansion ([5]) converges for small and 
moderate dissipation, up to e > 0.6, but breaks down 
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FIG. 1: Normal restitution coefficient e, as it follows from 
from Eq. for different values of A''^. Here m = 1 and 
(|ci2 • n|) = 4/v'27r. The latter value is the average dimen- 
sionless relative normal impact velocity with the averaging 
performed over a Maxwellian distribution. 



otherwise [3J|. If the dissipation is not large, one can 
consider only the first two non-trivial coefficients and 
03, which have been studied both analytically [3f. 33„,33| 
and numerically [H, ^ [11 . 

The theory developed for the case of constant e pre- 
dicts a rapid relaxation of to steady state values, which 
depend on e [H H^. The solution of the Boltzmann 
equation by the Direct Simulation Monte Carlo (DSMC) 
method demonstrated a good agreement between theo- 
retical and numerical results j34| . The analysis of the 
granular gas dynamics by the DSMC has however a seri- 
ous drawback - it does not take into account correlations 
between particles velocities which arise during system's 
evolution. It is not known a 'priori whether these cor- 
relations have an impact on the studied gas properties. 
Indeed, Molecular dynamics (MD) simulations which are 
free from the above shortcoming show at certain den- 
sities noticeable deviations of the numerically obtained 
coefficients a^ii) from the corresponding theoretical val- 
ues [1^, [11] . This discrepancy arises due to formation of 
spatial structures in a gas, which spontaneously appear 
when the gas evolves in a HCS. 
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Evolution of the velocity distribution function and re- 
spectively the Sonine coefficients is qualitatively different 
in a gas of viscoelastic particles. In contrast to the case of 
a constant e, where rapid relaxation is observed, the evo- 
lution of the Sonine coefficients 02 and 03 demonstrates 
a complicated non- monotonous time dependence jSTl [ssl ] . 
Typically, the magnitudes of the Sonine coefficients first 
increase with time, reach maximal values, and finally 
tend to zero, so that asymptotically the Maxwellian dis- 
tribution is recovered jSTl |38| . This simply follows from 
the fact that the collisions become more and more elastic 
as the average particles velocity decreases [cf. Eq. ([5])]. A 
theory describing evolution of 02 (t) and 03 (t) in a force- 
free granular gas of viscoelastic particles has been devel- 
oped in (ssj . however a corresponding theory for heated 
gases is still lacking. Moreover, the behavior of the ve- 
locity distribution function in a gas of viscoelastic parti- 
cles has never been studied by means of MD neither for 
force-free nor for driven systems. Hence, it is still not 
known how accurate is the theory for the velocity distri- 
bution function in this case and for which dissipation 6 
the Sonine series keeps its convergence. Another interest- 
ing question is the adequacy of the simplified model of a 
constant restitution coefficient for uniformly heated gases 
of viscoelastic particles: While it is known that model of 
constant e fails to describe a complicated evolution of a 
gas in a homogeneous cooling state [1^, it is still not 
known how good this model for driven systems. 



II. KINETIC THEORY OF UNIFORMLY 
HEATED GRANULAR GAS OF VISCOELASTIC 
PARTICLES 

We address a granular gas of identical viscoelastic par- 
ticles in a white- noise thermostat, that is, the gas under a 
uniform heating. The equation of motion for i-th particle 
then reads: 
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Here the force actin g be tween colliding i-th and j-th par- 
ticles is of the form [23] 

Fr" = E f/^cf + UpcII'^ e(o,)e,, (11) 



where p = ^niK (2/3)^^^ and Cij = i^ij — cr is the defor- 
mation of particles in a pairwise collision. It is expressed 



in terms of the inter-particle distance, 



and 



particle's diameter a. The unit vector = Tij/rij gives 
the direction of the inter-particle force and the Heaviside 
step-function 9(.t) indicates that the force acts only at a 
collision, when particles' centers approach closer than the 
diameter a. In Eq. (|lll) the summation over all grains is 
performed. For the external force F°''* we use the model 
of a Gaussian white noise: 



In the present study we perform thorough large-scale 
MD simulation studies of the velocity distribution func- 
tion in a gas of viscoelastic particles, both in a homoge- 
neous cooling state and under a uniform heating; we also 
develop a theory for /(v,t) for driven gases. We observe 
that our MD results for a2{t) and a^lt) are in excellent 
agreement with the analytical predictions of Ref. [s^ for 
the HCS, as well as with the results of the novel theory 
for driven systems, reported here. We present also simu- 
lation results for the Sonine coefficients 04 and 05 , which 
may be hardly obtained from the kinetic theory. The 
analysis of the series of these coefficients, 02, 03, 04, 05, 
allows to asses the convergence of the Sonine-polynomial 
expansion. The problem of an adequacy of the model 
of the effective constant restitution coefficient for driven 
gases of viscoelastic particles is also addressed. 



{Frit)) = 0, iK^m^Sit')) = m'mjS^pSit - 1') , 

(12) 

where a, /3 = x,y,z, and characterizes the magnitude 
of the stochastic force. 

We use event-driven simulations [s^, that is, we do 
not explicitly solve an equation of motion for colliding 
particles with the forces F'^°'\ but directly use the after- 
collision velocities, expressed in terms of the pre-coUision 
ones and the coefficient of restitution e [1, [s^ . If we fur- 
ther apply the hypothesis of molecular chaos, evolution 
of the velocity distribution function may be described us- 
ing the Boltzmann-Enskog (BE) equation, supplemented 
by the diffusive Fokker-Planck term, which corresponds 
to a uniform heating (2l| : 
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This paper is organized as follows. In the next Sec. II 
we develop a theory for evolution of the velocity distri- 
bution function in a heated gas of viscoelastic particles, 
that is, we obtain the time dependence of the Sonine co- 
efficients; the previous results for 02 (t) and a3(i) for a 
force-free gas are directly reproduced in a limiting case 
of vanishing heating. In Sec. Ill we report our MD re- 
sults and compare simulation data with the theoretical 
predictions. Finally, in Sec. IV we summarize our find- 
ings. 



Here /(/,/) is the collision integral and (72(0') is con- 
tact value of the pair correlation function that takes into 
account the excluded volume effect Q. For the scaling 
distribution function /(c, t) [cf. Eq. (jH])] , the above equa- 
tion reads 
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Here i{j,f^ — a '^v'^I{f,f) is the dimensionless 
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collision integral, defined as 



HfJ)^ J dc2 J dee{-c,2-e)\-c,2-e\ (15) 
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Here c" and Cj' are the (reduced) pre-coUision veloci- 
ties in the so-called inverse collision, which results in the 
(reduced) post-collision velocities Ci and C2 [1] . 

Now we substitute the Sonine-polynomial expansion 
for f{c,t) [Eq. ^] into Eq. and neglect all terms 
with p > 3. Multiplying the resulting equation with c\, 
c\ and , and integrating over Ci, we obtain 



dt 

d{c^) 
dt 

dt 



Here, 
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is the p-th moment of the dimensionless collision inte- 
gral. The moments fj,p can be calculated analytically up 
to 0{S^^) using a formula manipulation program as ex- 
plained in detail in Ref. [l]. They can be represented in 
the following form: 
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In the above equation we use the dimensionless time r = 
t/Tc{0), expressed in terms of the initial mean collision 
time Tc(0), where 



-\t) ^ 4V^g2{<j)<J^n 



T{t) 
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Eqs. ((20)) describe evolution of temperature and the ve- 
locity distribution function for a heated gas of viscoelas- 
tic particles. For K = 0, that is for = 0' Eqs. (gOl) 
reduce to corresponding equations of Ref. [s^ , which de- 
scribe evolution of 02 and 03 in a force-free gas. Note that 
Eqs. (PO)) also describe evolution of Sonine coefficients for 
a granular gas with a constant restitution coefficient. In 
this case the respective values of /ip for a constant e are 



to be used Ref. 34 



When a heated gas settles in a steady state, the gran- 
ular temperature and the velocity distribution function 
become time-independent. Then Eqs. pop yield the al- 
gebraic equations for the steady-state temperature and 
moments fip, p = 2, 4, 6: 



20 



k=0 



The numerical values of the coefficients M^'^l for fip with 
p = 2, 4, 6 arc given in Appendix. 

The Sonine coefficients are related to the moments of 
the reduced velocity as follows: 
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(1 + 02), 
(1 + 3a2 - aa) 



Hence using Eqs. one can obtain the following sys- 
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(23) 
(24) 
(25) 



Solving this system in the linear approximation with re- 
spect to the dissipation parameter S, one obtains the 
steady state Sonine coefficients 



-A26, 



A, 



500 



Ci ~ 0.435, 
(26) 

Ci ~ 0.078. 



(27) 

The above expressions for 02 and 03, together with 
Eq. (|19|) for /X2 yield the corresponding analytical ex- 
pression for the steady state temperature Tg.s. • 
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III. MOLECULAR DYNAMICS SIMULATIONS 
A. Force-free Gas 

Wc perform MD simulations [4l| , using an event-driven 
algorithm [s^ to investigate a system of inelastic gran- 
ular particles with the restitution coefficient ([5|), as it 
follows for viscoelastic particles. The grains are modeled 
as identical hard spheres of mass m = 1 and diameter 
a — 1. In event-driven simulations particles move freely 
between pairwise collisions, where the velocities of the 
particles are updated according to Eq. ([2]). Our system 
consists of = 2.048 x 10^ particles placed in a three- 
dimensional cubic box with the edge L = 418, which cor- 
responds to the number density of n = 0.028. Periodic 
boundary conditions are applied in all three directions. 

The system was initialized by placing randomly parti- 
cles in the box. For initial particles' velocities randomly 
directed vectors of the same length have been assigned. 
The initial total momentum mvi(O) was in this case 
zero with a high accuracy. We start with e = 1, which 
corresponds to elastic collisions, and allow the system to 
relax to a Maxwellian velocity distribution, which used 
to happen in a few collisions per particle. Hence, a ho- 
mogeneous system with a Maxwellian distribution was 
an initial condition for each MD run; we then evolve the 
system with a velocity-dependent restitution coefficient 
([5]). In simulations we used the initial granular temper- 
ature r(0) = 400/3 with the corresponding initial mean 
collision time Tc{Q) ~ 0.436. We have confirmed that the 
system remains approximately homogeneous during our 
simulations. All statistical quantities presented here have 
been obtained as averages over 25 independent runs. 

In Fig. [21 we present simulation results for the depen- 
dence of the reduced granular temperature u{t) on the 
dimensionless time r. In a HCS temperature decays ac- 
cording to the power-law m(t) - t^^/^ g[3j[3|]. The 
agreement between MD results (denoted by symbols) and 
theory (denoted by lines) is very good for both small and 
moderate values of the dissipative parameter 6; as ex- 
pected, the larger 6, the faster the gas cools. In Fig. [21 
we also show the results for the heated gas, which are 
discussed in the following subsection. Next, we analyze 
evolution of the velocity distribution, characterized by 
the time-dependent Sonine coefficients [cf. Eq. ([5])]. In 
Fig. [21 we show evolution of 02 and 03 for a force-free gas 
with (5 = 0.15. In order to calculate 02 and 03 numeri- 
cally, we compute the moments of the velocity distribu- 
tion function (c^*-') and use Eqs. (|19p . At initial time, 
the Soninc coefficients equal to zero: 02(0) — a^{0) — 0, 
which corresponds to a Maxwellian distribution. As it 
follows from Fig. [31 the absolute values of 02 and a^, ini- 
tially increase, reach a maximum, and eventually decay 
to zero. As the thermal velocity decreases during the 
gas cooling, the behavior of the system tends to that of 
a gas of elastic particles, that is, the velocity distribu- 
tion tends to a Maxwellian. The maximum values of \a2\ 
and |a3| increase with S, saturate at (5 ~ 0.15, and do 




FIG. 2: Evolution of the reduced granular temperature it(r) = 
r(r)/T(0) for the force-free and heated granular gases. The 
symbols denote simulation results, the lines - theoretical pre- 
dictions [the numerical solution of Eqs. (|20|l ]. Temperature 
of the force-free gas decreases asymptotically, r ^ 1, as 
u ~ r^^''^, while that of the heated gas first decreases and 
then saturates at a constant, steady state value Us.s.- Gray 
lines correspond to evolution of granular gas with the effective 
constant restitution coefficient Es.a. calculated from Eq. ^ 
with the same 5 and steady-state reduced temperature u^.a. 
as in the gas of viscoelastic particles, and with the averaged 
normal impact velocity, |ci2 • ej = 4/\/27r. 



not increase further; the location of these maxima shifts 
to later times with increasing S [38j . To investigate the 
convergence of the Sonine-polynomial expansion (|5]) we 
analyze the high-order Sonine coefficients, 04 and 05. In 
Fig. 31 the simulation data for the time dependence of 
0-2, 0.3, 0.4, as along with the analytical predictions for 02 
and 03 are shown for S = 0.1. Note, that 04 and 05 evolve 
in a similar fashion as 02 and 03, that is, their magni- 
tudes first increase, reach maximum and relax back to 
zero. Certainly, this is not surprising, since it is a simple 
manifestation of the fact that the system starts from the 
Maxwellian velocity distribution due to initial conditions 
and returns to it later, when the collisions become elas- 
tic. Again we see an excellent agreement between the 
simulation results and the theory for a2 and 03 . Further- 
more, 04 is an order of magnitude smaller than 02, while 
as is ~ 10^ times smaller than 02- This confirms the 
convergence of the Sonine-polynomial expansion for the 
velocity distribution, at least for these values of 6. 

In Fig. [HI we show the dependence of Sonince co- 
efficients on (dimensionless) time for larger dissipation, 
S — 0.3. The behavior of the Soninc coefficients is quali- 
tatively the same as for smaller 5. Moreover, the subse- 
quent coefficients ak decrease by an order of magnitude 
with increasing k, which again indicates the convergence 
of the Sonine-polynomial expansion for this dissipation. 
As the same time a significant discrepancy between the 
kinetic theory and the MD results is obvious. Fig. \5\ In 
the earlier work [s^, two of us have reported that the 
expansion ^ for e with Nk — 20 does not yield an ac- 
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FIG. 3: Evolution of a2 and 03 for the force-free and heated 
granular gases with 5 = 0.15. The symbols correspond to 
simulation results, the lines to our theory [numerical solution 
of Eqs. ([201)]. 



FIG. 4: Evolution of Sonine coefficients (k = 2, 3, 4, 5) for 
a force-free gas. Here the dissipation parameter 5 = 0.1. The 
symbols correspond to the MD results and the lines to the 
theory. 



ceptable accuracy for the Sonine coefficient 02 and as for 
5 > 0.3. The discrepancy seen in Fig. [S] between the 
theoretical curves and the simulation data supports this 
theoretical conclusion. 



B. Uniformly heated gas 

We study a uniformly-heated granular gas of viscoelas- 
tic particles using the algorithm suggested in Ref. [1^ 
(see also [H, During event-driven simulations, we 

heat the system after a time-step dt by adding to the 
velocity of each i-th particle a random increment, which 
mimics the heating by noise, 

Vi^a{t + dt) = v.^ait) + VrVdi^p , (28) 

where a = x,y, z. In Eq. (|28p. is a random number 
uniformly distributed within interval [^0.5, 0.5], and r is 
the amplitude of the noise, r = 12^q. After the change in 
velocities, the system is transferred to the centre-of-mass 



frame, 

1 ^ 

i=l 

which ensures conservation of the momentum. In our 
simulations dt = 0.1 and r = 0.1 were taken. We used 
the random kicking interval dt small, as compared to the 
mean free time. The number of particles and the sys- 
tem size were the same as in the force-free simulations; 
the initial conditions have been also prepared as in the 
force-free case. All statistical quantities were obtained as 
averages over 10 independent runs. The results are shown 
in Figs. El El [7] and El As it follows from Fig. El during 
the early stages of cooling, the temperature of a heated 
granular gas evolves in the same way as the respective 
force-free gas with the same 6. However, while the tem- 
perature of the heated gas settles at a steady-state value, 
the force-free gas continues to cool down. The results of 
our kinetic theory for u{t) perfectly agree with the sim- 
ulation data. In Fig. |31 we plot 02, 03 as functions of 
the dimcnsionless time r for both force-free and heated 
gas with 5 = 0.15. For the initial stage of evolution, at 
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FIG. 5: Evolution of a2, as, 04, as for a free granular gas with 
S = 0.3. 



T < 100, the Sonine coefScients for heated and force-free 
gases almost coincide. At later times, the coefficients of 
the heated gas settle at the steady-state [cf. Eqs. (PS)) - 
(P?])]. whereas these for the freely evolving gas continue 
to decay. Again, as for the force-free gases, our theory is 
again in excellent agreement with computer simulations. 

Figure [B] demonstrates the convergence of the Sonine- 
polynomial expansion for heated granular gases. As in 
Fig. m we plot here the dependence of 02, 03, 04 and 05 
on r for S = 0.1. We see that successive Sonine coeffi- 
cients decrease by an order of magnitude, which clearly 
indicates the convergence of the Sonine series. Finally, in 
Fig. [71 we plot 02, 03, 04 and 05 for the heated gas with 
relatively large dissipation, S = 0.3. In this case one can 
see a noticeable discrepancy between our theory and the 
numerical results (cf. Fig. [5] for the freely- evolving gas). 
Nevertheless, smallness of the high-order coefficients 04 
and 05 seemingly manifests the convergence of the So- 
nine series. Interestingly, although the simulation data 
deviate from the theoretical curves during the relaxation 
to the steady state, the numerical and theoretical values 
of the Sonine coefficients 02 and as in the steady state 
are very close, see Fig. [7] for r > 10"^. 



FIG. 6: Evolution of Sonine coefficients at {k = 2, 3, 4, 5) for 
a uniformly-heated gas. Here, the dissipation parameter 5 = 
0.1. Gray lines correspond to evolution of granular gas with 
effective constant restitution coeffcient, calculated according 
to Eq. ((si with same 5 and steady-state relative temperature 
u, as in the viscoelastic case, and {|ci2 ■ n|) = 4/v'27r. 



C. Comparison with the case of a constant 
restitution coefficient 

As we have already mentioned in the Introduction, the 
behavior of a granular gas with a constant restitution 
coefficient in a homogeneous cooling state qualitatively 
differs from that of a gas of viscoelastic particles. The 
reason for this is the permanent cooling of the gas, which 
causes the permanent increase of the effective restitution 
coefficient. Contrary, uniformly heated gases evolve to a 
steady state, with constant temperature Tg.s. where an ef- 
fective restitution coefficient is kept constant. Therefore, 
it is natural to assume a somewhat close behavior of the 
two gases, with the same initial temperature Tq, number 
density and heating rate, provided the effective restitu- 
tion coefficient Eg.s. is appropriately chosen. To chose Eeff 
we assume that it should be equal to the restitution coef- 
ficient of viscoelastic particles colliding with the relative 
velocity, equal thermal velocity at the steady state Tg.s., 
with the average normal component of the impact veloc- 
ity (|ci2 • n|). That is, to compute e^s we exploit Eq. ^ 
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FIG. 7: Evolution of 02, as, a4, as for a uniformly-heated gas 
with S = 0.3. Gray lines correspond to evolution of granular 
gas with effective constant restitution coeffcient, calculated 
according to Eq. ([SJ with same 5 and steady-state relative 
temperature u, as in the viscoelastic case, and (|ci2 ■ n|) — 

with u{t) = Us.s. = Ts.s./Tf) and (|ci2 ■ n|) = 4/-\/27r, with 
the corresponding dissipation 6. Fig. [5] demonstrates al- 
most coinciding behavior of the granular temperatures 
in the two gases. Evolution of temperature towards the 
steady-state in the in the gas with the constant restitu- 
tion coefficient e = Coff is practically indistinguishable 
from that of the gas of viscoelastic particles; a small dis- 
crepancy may be seen only for intermediate times. 

Nevertheless, a quantitative agreement in the evolu- 
tion of the velocity distribution function of the two sys- 
tems is lacking; only qualitative agreement is observed 
in the time dependence of the Sonine coefficients. Still 
however, when the dissipation is not large, S = 0.1, the 
steady state values of 02 and 03 are almost the same for 
the two gases, Fig. [HI For larger dissipation, S = 0.3, 
not only evolution of these coefficients is rather differ- 
ent, but also their steady-state limits, Fig. [71 Hence we 
conclude that while for small dissipation it is possible to 
use an effective constant restitution coefficient to model 
simultaneously temperature and the form of the veloc- 
ity distribution function (quantified by 02 and 03) in a 



steady state, this can not be done for large dissipations. 

Therefore, as it follows from our study a simplified col- 
lision model with an effective constant restitution coeffi- 
cient may be used only to describe some averaged, that 
is, macroscopic characteristics, like temperature, but not 
microscopic, as the velocity distribution function. 

IV. CONCLUSION 

To date, most studies of granular gases have been fo- 
cused on systems with a constant restitution coefficient 
e, although the basic physics requires impact- velocity de- 
pendence of e. Here we have studied numerically, by 
means of Molecular Dynamics, and theoretically a gran- 
ular gas of particles with the impact- velocity dependent 
restitution coefficient, as it follows from the simplest first- 
principle model of viscoelastic spheres. To explore evolu- 
tion of granular temperature as well as the velocity distri- 
bution function we performed large-scale MD simulations 
of a force-free gas and of a gas under a uniform heating. 
We also developed a theory of the velocity distribution 
function for a driven gas of viscoelastic particles. The ve- 
locity distribution function deviates from a Maxwellian 
distribution and deviations are quantified by Sonine coef- 
ficients in the Sonine polynomials expansion. In our MD 
simulations we have found a few first Sonine coefficients 
02, 03, a4, as and obtained theoretical predictions for 02 
and a^. We noticed that each next order Sonine coeffi- 
cient decreases by an order of magnitude with respect to 
the previous one, that is, 03 is approximately 10 times 
smaller than 02, while 04 is 10 times smaller than 03 but 
10 time larger than 05. This implies the convergence of 
the Sonine series for both heated and force-free gases in 
the range of dissipation studied. 

We observed that predictions of our theory are in a 
very good agreement with the MD results for granular 
temperature for both small and large dissipation, quan- 
tified by the dissipative parameter i5. We also found that 
if dissipation is not large, (S = 0.1), the dependence on 
time of the Sonine coefficient, obtained in simulations is 
in excellent agreement with the theoretical predictions 
for both cases of a homogeneous cooling state (existed 
theory) and of a uniform heating (the novel theory); at 
the same time for large dissipation {S = 0.3) a notice- 
able deviation between the theory and MD results are 
observed for. 

We also analyze the accuracy of the simplified model of 
the effective constant restitution coefficient Scs for uni- 
formly heated gases. We reveal, that while this model 
may be used to describe the temperature of a gas, it fails 
to describe quantitatively both temperature and the ve- 
locity distribution function simultaneously, especially for 
large dissipations. 
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V. APPENDIX 

Here we present numerical values of the coefficients M^'^^ for p-th moments of the dimensionless collision integral 

fip = — J dcd'I{f, /), whereo = 2, 4, 6 (Eq. If 9p . They can be calculated using a formula manipulation program as 
explained in detail in Rcf. 



TABLE I: Numerical coefficients for fi2- 



k 








^4 22 


^"2. 32 


"-1-2, 23 























1 




















2 


6.485666782 


1.55656004 


0.103770674 


0.05447964 


0.012347568 


0.043583201 


3 




















4 


-9.285794108 


-2.762523735 


-0.138126145 


-0.067336539 


-0.014198354 


-0.051624196 


5 


-1.801263542 


-0.591039599 


-0.024626646 


-0.011543787 


-0.002345029 


-0.008657705 


6 


10.39566154 


3.742438128 


0.124747864 


0.0561366 


0.010977733 


0.041168185 


7 


5.911161093 


2.32382521 


0.058095712 


0.025053679 


0.004712918 


0.017955684 


8 


-10.43948573 


-4.462880197 


-0.074381459 


-0.030682122 


-0.005549036 


-0.02147679 


9 


-10.53072647 


-4.877042711 


-0.040642187 


-0.016002713 


-0.002778439 


-0.01093409 


10 


8.770061385 


4.385030729 


-6.28319E-08 


1.5708E-07 


-6.9115E-07 


-1.09956E-06 


11 


14.12940508 


7.603386089 


-0.063361675 


-0.022572577 


-0.003606668 


-0.014669549 


12 


-4.537693883 


-2.620518195 


0.043675344 


0.014740329 


0.00225554 


0.009334873 


13 


-16.3720861 


-10.11999574 


0.25299992 


0.080642707 


0.011807124 


0.04972868 


14 


-1.57951827 


-1.042482063 


0.034749403 


0.010424777 


0.001458779 


0.006254857 


15 


16.74058856 


11.77072632 


-0.490447192 


-0.137938722 


-0.018439173 


-0.080466684 


16 


8.091094766 


6.048093381 


-0.302404622 


-0.079381265 


-0.010120669 


-0.044983133 


17 


-14.28935468 


-11.33324447 


0.661105967 


0.161144633 


0.019577921 


0.088628216 


18 


-13.91371204 


-11.68751806 


0.779168081 


0.175313108 


0.020253112 


0.093501317 


19 


8.817864446 


7.831365892 


-0.587352513 


-0.12114108 


-0.013292359 


-0.062590324 


20 


18.12640747 


16.99350704 


-1.416125536 


-0.265523424 


-0.027660643 


-0.132762689 



TABLE II: Numerical coefficients for /.t4. 



k 






^"4. 3 


^"4, 22 


Mn2 


-'"4, 23 








10.0265131 


-2.506628274 


0.313328534 


0.048957583 


0.156664267 


1 




















2 


36.31973398 


46.85245668 


-6.501232013 


-0.294188929 


-0.014744639 


-0.137283095 


3 




















4 


-71.50061454 


-100.6607515 


16.08697411 


0.588254117 


0.028578585 


0.254896238 


5 


-10.35726536 


-15.38813811 


2.65528065 


0.077425377 


0.003372684 


0.032260122 


6 


116.2110953 


170.7289911 


-29.24466734 


-0.978197119 


-0.053925433 


-0.391290111 


7 


45.69327521 


73.08622368 


-13.9966349 


-0.27504694 


-0.011868725 


-0.105434433 


8 


-169.5076903 


-260.3837618 


47.12363979 


1.010942243 


0.055144437 


0.370717555 


9 


-117.6984195 


-196.0891402 


39.54174543 


0.343236574 


0.015776747 


0.120152601 


10 


219.5627451 


354.2742165 


-67.35573679 


9.42478E-07 


-9.59128E-05 


2.19911E-06 


11 


234.6065359 


406.4486348 


-85.0337586 


0.942565153 


0.041890374 


0.298516588 


12 


-240.6229512 


-406.3774129 


79.91582806 


-2.489342054 


-0.113668109 


-0.746801174 


13 


-398.7161877 


-724.5838699 


157.8564548 


-6.21062323 


-0.248019053 


-1.759667007 


14 


198.2681578 


345.5995343 


-67.77893928 


5.334853916 


0.215424796 


1.422686368 


15 


597.9839776 


1147.255623 


-261.5391558 


19.23940936 


0.663323094 


4.809946224 


16 


-57.48669811 


-83.77829734 


7.187455294 


-3.892354379 


-0.14905002 


-0.908182286 


17 


-797.173725 


-1618.265223 


386.3071566 


-43.41915585 


-1.250041435 


-9.407465346 


18 


-208.1541912 


-471.7945326 


132.0933869 


-11.561012 


-0.264278169 


-2.312272525 


19 


935.7389044 


2010.206189 


-501.6884258 


77.95690697 


1.814066361 


14.29254849 


20 


604.5584047 


1387.679757 


-379.5240628 


55.87757953 


1.100691783 


9.313324364 
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TABLE III: Numerical coefficients for /ig. 



k 


», r(*-) 
^"6, 




^"6, 3 




-"^6, 32 


-"^6, 23 








112.7982724 


-84.5987043 


-4.934924416 


0.022030913 


-0.822487403 


1 




















2 


209.9410334 


633.7845455 


-245.0212299 


16.75084546 


0.263561766 


2.233405168 


3 




















4 


-525.042013 


-1718.358426 


717.5389515 


-48.5267813 


-0.447724986 


-4.852819306 


5 


-61.2148156 


-203.8361017 


86.29524105 


-7.201298665 


-0.061076884 


-0.600095271 


6 


1097.056804 


3731.587123 


-1607.525698 


133.8829794 


0.639814941 


8.92527426 


7 


338.2152097 


1206.681042 


-544.0572385 


46.03791709 


0.181904962 


2.302036496 


8 


-2042.435348 


-7140.741695 


3126.21351 


-347.473825 


-0.750086425 


-11.58235962 


9 


-1116.580388 


-4130.985458 


1913.431002 


-196.4184963 


-0.205720952 


-3.273691172 


10 


3411.978481 


12219.71999 


-5395.763075 


806.6588655 


-0.002230531 


-0.001201659 


11 


2856.108247 


10832.88409 


-5076.69307 


683.3582205 


-0.646810061 


-11.39021812 


12 


-5055.198635 


-18485.31829 


8165.55863 


-1634.92685 


3.291775252 


54.49778445 


13 


-6187.93005 


-24012.30452 


11320.12826 


-2006.585914 


5.40301205 


100.3307334 


14 


6435.915735 


23886.43086 


-10416.97776 


2818.339574 


-10.83459224 


-187.8895259 


15 


11807.48498 


46936.94502 


-22224.52202 


5062.759585 


-21.37481894 


-421.907027 


16 


-6433.236775 


-23805.35645 


9879.364245 


-3867.079426 


21.1398954 


386.7062636 


17 


-20203.41898 


-82418.70275 


39185.83747 


-11160.63025 


60.2704598 


1302.074732 


18 


3203.277562 


10083.97492 


-2516.563404 


3251.164745 


-23.50300074 


-433.4876022 


19 


31188.58356 


130743.2016 


-62400.7815 


21721.86788 


-133.1393616 


-3258.298424 


20 


5784.30056 


29433.94954 


-17812.70067 


2307.910281 


-6.449835295 


-384.6515338 
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